# Author: Patrick Cunha Silva
# Project: Janusz, Cunha, and Junqueira
# Goal: SI for General Tendencies
# First Version: Dec 20, 2023
# This Version: Dec 20, 2023

rm(list = ls())

# Load data
load(file = here::here("Data", "descrData.RData"))

# To count the number of observations
descrData$x <- 1


#######################################
#######################################
############### Women #################
#######################################
#######################################

#### Raw Number
# Calculate Statistics
mean_women_ver <- with(descrData, aggregate(woman_ver, list(ANO_ELEICAO, woman_pref), mean))
sd_women_ver <- with(descrData, aggregate(woman_ver, list(ANO_ELEICAO, woman_pref), sd))
n_women_ver <- with(descrData, aggregate(x, list(ANO_ELEICAO, woman_pref), sum))
# Calculate SE 
se <- sd_women_ver$x/sqrt(n_women_ver$x)
# Add CI
mean_women_ver$ub <- mean_women_ver$x + qnorm(0.975) * se
mean_women_ver$lb <- mean_women_ver$x + qnorm(0.025) * se
# Sort data
mean_women_ver <- mean_women_ver[order(mean_women_ver$Group.1, mean_women_ver$Group.2),]

##### Ratio
# Calculate Statistics
mean_women_ratio <- with(descrData, aggregate(woman_ver_ratio, list(ANO_ELEICAO, woman_pref), mean))
sd_women_ratio <- with(descrData, aggregate(woman_ver_ratio, list(ANO_ELEICAO, woman_pref), sd))

# Calculate SE 
se <- sd_women_ratio$x/sqrt(n_women_ver$x)
# Add CI
mean_women_ratio$ub <- mean_women_ratio$x + qnorm(0.975) * se
mean_women_ratio$lb <- mean_women_ratio$x + qnorm(0.025) * se
# Sort data
mean_women_ratio <- mean_women_ratio[order(mean_women_ratio$Group.1, mean_women_ratio$Group.2),]

##### %
# Calculate Statistics
mean_women_pct <- with(descrData, aggregate(woman_ver_pct, list(ANO_ELEICAO, woman_pref), mean))
sd_women_pct <- with(descrData, aggregate(woman_ver_pct, list(ANO_ELEICAO, woman_pref), sd))

# Calculate SE 
se <- sd_women_pct$x/sqrt(n_women_ver$x)
# Add CI
mean_women_pct$ub <- mean_women_pct$x + qnorm(0.975) * se
mean_women_pct$lb <- mean_women_pct$x + qnorm(0.025) * se
# Sort data
mean_women_pct <- mean_women_pct[order(mean_women_pct$Group.1, mean_women_pct$Group.2),]


#######################
######Create Plots#####
#######################

pdf(here::here('plots_tables', 'WomenDescriptive.pdf'), height = 8, width = 8)
# Create plot
par(mar = c(5, 6, 2, 1))
myP <- barplot(mean_women_ver$x, col = c('grey80', 'grey20'),
               ylim = c(0, 45),
               ylab = 'Average Ratio of Women Running', 
               cex.lab = 2, cex.axis = 2, 
               border = FALSE)
indexer1 <- seq(1, 10, 2)
indexer2 <- seq(2, 10, 2)
atP <- c()
for(i in 1:5) atP <- c(atP, (myP[indexer1[i], 1] + myP[indexer2[i], 1])/2)
axis(1, at = atP, labels = unique(mean_women_ver$Group.1), cex.axis = 2, tick = F)
arrows(myP, mean_women_ver$ub, myP, mean_women_ver$lb, angle = 90, code = 3 , length=0.1)
# legend(x = atP[4], y = 45, legend = c('Man Mayor', 'Woman Mayor'),
#        fill = c('grey85', 'grey25'), cex = 1.75, bty = 'n',
#        border = FALSE)

# Create plot
par(mar = c(5, 6, 2, 1))
myP <- barplot(mean_women_ratio$x, col = c('grey80', 'grey20'),
               ylim = c(0, 3.5),
               ylab = 'Average Ratio of Women Running', 
               cex.lab = 2, cex.axis = 2, 
               border = FALSE)
indexer1 <- seq(1, 10, 2)
indexer2 <- seq(2, 10, 2)
atP <- c()
for(i in 1:5) atP <- c(atP, (myP[indexer1[i], 1] + myP[indexer2[i], 1])/2)
axis(1, at = atP, labels = unique(mean_women_ratio$Group.1), cex.axis = 2, tick = F)
arrows(myP, mean_women_ratio$ub, myP, mean_women_ratio$lb, angle = 90, code = 3 , length=0.1)
# legend(x = atP[4], y = 3.5, legend = c('Man Mayor', 'Woman Mayor'),
#        fill = c('grey85', 'grey25'), cex = 1.75, bty = 'n',
#        border = FALSE)


# Create plot
par(mar = c(5, 6, 2, 1))
myP <- barplot(mean_women_pct$x, col = c('grey80', 'grey20'),
               ylim = c(0, 40),
               ylab = 'Average % of Women Running', 
               cex.lab = 2, cex.axis = 2, 
               border = FALSE)
indexer1 <- seq(1, 10, 2)
indexer2 <- seq(2, 10, 2)
atP <- c()
for(i in 1:5) atP <- c(atP, (myP[indexer1[i], 1] + myP[indexer2[i], 1])/2)
axis(1, at = atP, labels = unique(mean_women_pct$Group.1), cex.axis = 2, tick = F)
arrows(myP, mean_women_pct$ub, myP, mean_women_pct$lb, angle = 90, code = 3 , length=0.1)
legend(x = myP[1]-0.5, y = 40, legend = c('Man Mayor', 'Woman Mayor'),
       fill = c('grey85', 'grey25'), cex = 1.75, bty = 'n',
       border = FALSE)
dev.off()
                  

#################################################
#################################################
############### Afro-Brazilians #################
#################################################
#################################################

descrData <- subset(descrData, ANO_ELEICAO %in% c(2016, 2020))
#### Raw Number
# Calculate Statistics
mean_afroBR_ver <- with(descrData, aggregate(afroBR_ver, list(ANO_ELEICAO, afroBR_pref), mean))
sd_afroBR_ver <- with(descrData, aggregate(afroBR_ver, list(ANO_ELEICAO, afroBR_pref), sd))
n_afroBR_ver <- with(descrData, aggregate(x, list(ANO_ELEICAO, afroBR_pref), sum))
# Calculate SE 
se <- sd_afroBR_ver$x/sqrt(n_afroBR_ver$x)
# Add CI
mean_afroBR_ver$ub <- mean_afroBR_ver$x + qnorm(0.975) * se
mean_afroBR_ver$lb <- mean_afroBR_ver$x + qnorm(0.025) * se
# Sort data
mean_afroBR_ver <- mean_afroBR_ver[order(mean_afroBR_ver$Group.1, mean_afroBR_ver$Group.2),]

##### Ratio
# Calculate Statistics
mean_afroBR_ratio <- with(descrData, aggregate(afroBR_ver_ratio, list(ANO_ELEICAO, afroBR_pref), mean))
sd_afroBR_ratio <- with(descrData, aggregate(afroBR_ver_ratio, list(ANO_ELEICAO, afroBR_pref), sd))

# Calculate SE 
se <- sd_afroBR_ratio$x/sqrt(n_afroBR_ver$x)
# Add CI
mean_afroBR_ratio$ub <- mean_afroBR_ratio$x + qnorm(0.975) * se
mean_afroBR_ratio$lb <- mean_afroBR_ratio$x + qnorm(0.025) * se
# Sort data
mean_afroBR_ratio <- mean_afroBR_ratio[order(mean_afroBR_ratio$Group.1, mean_afroBR_ratio$Group.2),]

##### %
# Calculate Statistics
mean_afroBR_pct <- with(descrData, aggregate(afroBR_ver_pct, list(ANO_ELEICAO, afroBR_pref), mean))
sd_afroBR_pct <- with(descrData, aggregate(afroBR_ver_pct, list(ANO_ELEICAO, afroBR_pref), sd))

# Calculate SE 
se <- sd_afroBR_pct$x/sqrt(n_afroBR_ver$x)
# Add CI
mean_afroBR_pct$ub <- mean_afroBR_pct$x + qnorm(0.975) * se
mean_afroBR_pct$lb <- mean_afroBR_pct$x + qnorm(0.025) * se
# Sort data
mean_afroBR_pct <- mean_afroBR_pct[order(mean_afroBR_pct$Group.1, mean_afroBR_pct$Group.2),]


#######################
######Create Plots#####
#######################

pdf(here::here('plots_tables', 'afroBRDescriptive.pdf'), height = 8, width = 8)
# Create plot
par(mar = c(5, 6, 2, 1))
myP <- barplot(mean_afroBR_ver$x, col = c('grey50', 'grey80', 'grey20'),
               ylim = c(0, 80),
               ylab = 'Average Ratio of Afro-Brazilians Running', 
               cex.lab = 2, cex.axis = 2, 
               border = FALSE)
axis(1, at = c(myP[1], sum(myP[-1])/2), labels = unique(mean_afroBR_ver$Group.1), cex.axis = 2, tick = F)
arrows(myP, mean_afroBR_ver$ub, myP, mean_afroBR_ver$lb, angle = 90, code = 3 , length=0.1)
# legend(x = myP[2], y = 80, legend = c('All Mayors', 'Non AB Mayor', 'AB Mayor'),
#        fill = c('grey50', 'grey80', 'grey20'), cex = 1.75, bty = 'n',
#        border = FALSE)

# Create plot
par(mar = c(5, 6, 2, 1))
myP <- barplot(mean_afroBR_ratio$x, col = c('grey50', 'grey80', 'grey20'),
               ylim = c(0, 6),
               ylab = 'Average Ratio of Afro-Brazilians Running', 
               cex.lab = 2, cex.axis = 2, 
               border = FALSE)
axis(1, at = c(myP[1], sum(myP[-1])/2), labels = unique(mean_afroBR_ratio$Group.1), cex.axis = 2, tick = F)
arrows(myP, mean_afroBR_ratio$ub, myP, mean_afroBR_ratio$lb, angle = 90, code = 3 , length=0.1)
# legend(x = myP[1]-0.5, y = 6, legend = c('All Mayors', 'Non AB Mayor', 'AB Mayor'),
#        fill = c('grey50', 'grey80', 'grey20'), cex = 1.75, bty = 'n',
#        border = FALSE)

# Create plot
par(mar = c(5, 6, 2, 1))
myP <- barplot(mean_afroBR_pct$x, col = c('grey50', 'grey80', 'grey20'),
               ylim = c(0, 80),
               ylab = 'Average % of Afro-Brazilians Running', 
               cex.lab = 2, cex.axis = 2, 
               border = FALSE)
axis(1, at = c(myP[1], sum(myP[-1])/2),  labels = unique(mean_afroBR_pct$Group.1), cex.axis = 2, tick = F)
arrows(myP, mean_afroBR_pct$ub, myP, mean_afroBR_pct$lb, angle = 90, code = 3 , length=0.1)
legend(x = myP[1]-0.5, y = 80, legend = c('All Mayors', 'Non AB Mayor', 'AB Mayor'),
       fill = c('grey50', 'grey80', 'grey20'), cex = 1.75, bty = 'n',
       border = FALSE)
dev.off()
